function h = lpcvd_1DSurface_ConstructFromSnapshot(filename,LatticeSize)
% DESCRIPTION:
%  Find surface profile from snapshot file for 1D triangular lattice

fid = fopen(filename);
data = textscan(fid,'%f %f','HeaderLines',3);
fclose(fid);
layer = data{1}; % NOTE: the value in layer starts from 0
site = data{2};  % NOTE: the value in site is from 0 to (LatticeSize-1)

h = ones(1,2*LatticeSize);
TopLayer = max(layer);
if mod(TopLayer,2)==0
    TopEvenLayer = TopLayer;
    TopOddLayer = TopLayer-1;
else
    TopEvenLayer = TopLayer-1;
    TopOddLayer = TopLayer;
end
h(1:2:2*LatticeSize) = TopEvenLayer*h(1:2:2*LatticeSize);
h(2:2:2*LatticeSize) = TopOddLayer*h(2:2:2*LatticeSize);
% Find the top particles
for i = length(layer):-1:1
    if mod(layer(i),2)==0
        indx = (site(i))*2+1;
    else
        indx = (site(i)+1)*2;
    end
    if(h(indx)==layer(i))
        h(indx)=h(indx)-2;
    end    
end
hts_half = h(1:2:2*LatticeSize);
hts = h(2:2:2*LatticeSize);
NEXT_ST = [2:LatticeSize,1];
PREV_ST = [LatticeSize,1:(LatticeSize-1)];
for i = 0:(LatticeSize-1)
    h(2*i+1)=hts_half(i+1);
    h(2*i+2)=hts(i+1);
    if((hts_half(i+1)>hts(i+1))&&(hts_half(NEXT_ST(i+1))>hts(i+1)))
        h(2*i+2)=(hts_half(i+1)+hts_half(NEXT_ST(i+1)))/2;
    elseif((hts(PREV_ST(i+1))>hts_half(i+1))&&(hts(i+1)>hts_half(i+1)))
        h(i*2+1)=(hts(PREV_ST(i+1))+hts(i+1))/2;
    end        
end 
x = 0:0.5:(LatticeSize-0.5);
% for i = 1:2*LatticeSize
%     rectangle('position',[x(i),h(i),1,1],'FaceColor',[1,0,0],'EdgeColor',[0,0,0]);
% end
plot(0:0.5:(LatticeSize-0.5),h+0.65,'k-','LineWidth',1.5);
% if h(1)<h(2) && h(1)<h(2*LatticeSize)
%     h(1) = (h(2)+h(2*LatticeSize))/2;
% end
% if h(end)<h(1) && h(end)<h(end-1)
%     h(end) = (h(end-1)+h(1))/2;
% end
% for i = 2:1:(2*LatticeSize-1)
%     if (h(i)<h(i+1)&&h(i)<h(i-1))
%         h(i)=(h(i+1)+h(i-1))/2;
%     end
% end